%% ��������
close all
load data10;
k=[0.91902,0.040072];
ri=356;
%% ͨ��
%% ���Ƚ���AR(1)ģ��
qxs=reshape(qx,96,ri);
cfs=reshape(cf,96,ri);
arbc(1,:)=es(end,1:ri-1)*k(1)+k(2);
for i=2:96
    arbc(i,:)=arbc(i-1,:)*k(1)+k(2);
end
qx_ar=qxs(:,2:end)+arbc;
qx_ar(qx_ar<0)=0;
error_ar=cfs(:,2:end)-qx_ar;
for i=1:96
    e_ar(i,:)=mae(error_ar(i,:));
end

%% �����������Իع�ģ��
inputSeries=num2cell(qx');
targetSeries=num2cell(e');

inputDelays = 0;%�ⲿ�����0���ӳ�
feedbackDelays = 1:1;%�Իع��1���ӳ�
hiddenLayerSize = 8;
net = narxnet(inputDelays,feedbackDelays,hiddenLayerSize);

%% ��������Ԥ������������,�����䲻��Ԥ����������,���Ǹ�һ���������
net.inputs{1}.processFcns = {'removeconstantrows','mapminmax'};
net.inputs{2}.processFcns = {'removeconstantrows','mapminmax'};

%% ʱ����������׼������
[inputs,inputStates,layerStates,targets] = preparets(net,inputSeries,{},targetSeries);%������inputs,��inputSeries,targetSeries����cell��ʽ��������2��,��֪����ô���ɵ�

%% ѵ�����ݡ���֤���ݡ��������ݻ���
net.divideFcn = 'dividerand';  
net.divideMode = 'value';  
net.divideParam.trainRatio = 70/100;
net.divideParam.valRatio = 15/100;
net.divideParam.testRatio = 15/100;

%% ����ѵ�������趨
net.trainFcn = 'trainlm';  % Levenberg-Marquardt

%% �����趨
net.performFcn = 'mse';  % Mean squared error

%% ��ͼ�����趨
net.plotFcns = {'plotperform','plottrainstate','plotresponse', ...
  'ploterrcorr', 'plotinerrcorr'};

%% ����ѵ��
[net,tr] = train(net,inputs,targets,inputStates,layerStates);

%% �������
outputs = net(inputs,inputStates,layerStates);
errors = gsubtract(targets,outputs);
performance = perform(net,targets,outputs)

%% ����ѵ��������֤�������Լ����
trainTargets = gmultiply(targets,tr.trainMask);
valTargets = gmultiply(targets,tr.valMask);
testTargets = gmultiply(targets,tr.testMask);
trainPerformance = perform(net,trainTargets,outputs)
valPerformance = perform(net,valTargets,outputs)
testPerformance = perform(net,testTargets,outputs)

%% ����ѵ��Ч�����ӻ�
% figure, plotperform(tr)
% figure, plottrainstate(tr)
% figure, plotregression(targets,outputs)
% figure, plotresponse(targets,outputs)
% figure, ploterrcorr(errors)
% figure, plotinerrcorr(inputs,errors)

%% close loopģʽ��ʵ��-�Ǿ�ֵ����
% ����NARX������ģʽ
narx_net_closed = closeloop(net);
% view(net)
% view(narx_net_closed)

for i=1:ri-1
    phInputs_c=num2cell(qx(96*i:96*(i+1))');
    PhTargets_c=num2cell(e(96*i:96*(i+1))');

    [p1,Pi1,Ai1,t1] = preparets(narx_net_closed,phInputs_c,{},PhTargets_c);
    % �������
    yp1 = narx_net_closed(p1,Pi1,Ai1);
    qx_narx(:,i)=qxs(:,i+1)+cell2mat(yp1)';
%     figure;
%     plot([cell2mat(yp1)' cell2mat(t1)'])
end
qx_narx(qx_narx<0)=0;
error_narx=cfs(:,2:end)-qx_narx;
for i=1:96
    e_narx(i,:)=mae(error_narx(i,:));
end


% save data10_deal
%% close loopģʽ��ʵ��-��ֵ����,��ѵ�����е�Ч��,���ǶԲ��Լ��Ͳ�һ������,���Ըо����ǲ��ܼ���
% es_narx=[es(end,1:ri-1);es(:,2:ri)];
% qx_narx=[qxs(end,1:ri-1);qxs(:,2:ri)];
% 
% for i=1:96
%     for j=1:ri-1
%         phInputs_c=num2cell(qx_narx(i:i+1,j)');
%         if i==1
%             PhTargets_c=num2cell(es_narx(1:2,j)');
%         else
%             PhTargets_c=[yp1(i-1,j),num2cell(es_narx(i+1,j))];
%         end
%        [p1,Pi1,Ai1,t1(i,j)] = preparets(narx_net_closed,phInputs_c,{},PhTargets_c);
%        yp1(i,j) = narx_net_closed(p1,Pi1,Ai1);
%     end
%     error_narx2(i,:)=cell2mat(yp1(i,:))'-cell2mat(t1(i,:))';
%     m_narx(i,:)=mean(error_narx2(i,:));
%     yp1(i,:)=num2cell(cell2mat(yp1(i,:))-m_narx(i,:));
%     error_narxm(i,:)=error_narx2(i,:)-m_narx(i,:);
% end
% 
% for i=1:96
%     e_narxm(i,:)=mae(error_narxm(i,:));
% end
%% ��������ķ���Ԥ��ֵ����ܶȽ��к��ܶȹ���,����������,��Ϊ��ͬ�״εĺ��ܶȹ��ƻ��ͬ,ȷ��������Ԥ��
qwb=cell(1,96);
for j=1:96
    enarx=error_narx(j,:)';
    qxnarx=qx_narx(j,:)';
    gridx1=min(enarx):.1:max(enarx);%�������������,����0.1
    gridx2=min(qxnarx):.1:max(qxnarx);%���ַ���������,����0.1
    [x1,x2]=meshgrid(gridx1, gridx2);
    x1=x1(:);
    x2=x2(:);
    xi=[x1 x2];%����4��,Ϊ�̶�ʽ,���γ�2ά����,ÿһ����һ�������
    X=[enarx,qxnarx];%�������������ٽ�������,�γ�ʵ�����ݵ�ռ�
    [gl,bl]=ksdensity(X,xi);
    if j==48 
        figure;
        ksdensity(X,xi);%��������ֱ�ӻ�ͼ,ʵ��ʵ�ֲ��˸����������ݻ�,��ɫ���Ǻڵ�
%         Rcolorbar;
        x1=xlabel('Wind Speed Error (m/s)');
        y1=ylabel('Wind Speed (m/s)');
        set(x1,'Rotation',13);%x��������ת
        set(y1,'Rotation',-16.5);%y��������ת
        zlabel('Kernel Density');
        set(gca, 'FontSize', 30);%�ֺ���Ϊ30
    end
    % �γɸý׵ķ�λ����
    qw=[];
    for i=min(qxnarx):0.1:max(qxnarx)
        clear he
        I=find(abs(bl(:,2)-i)<1e-2);
        zh=sum(gl(I));
        he(1)=gl(I(1));
        for k=2:length(I)
            he(k,:)=he(k-1)+gl(I(k));
        end
        s90=bl(I(find(he>0.95*zh,1)),1);
        x90=bl(I(find(he>0.05*zh,1)),1);
        s80=bl(I(find(he>0.9*zh,1)),1);
        x80=bl(I(find(he>0.1*zh,1)),1);
        s70=bl(I(find(he>0.85*zh,1)),1);
        x70=bl(I(find(he>0.15*zh,1)),1);
        qw=[qw;i,s90,x90,s80,x80,s70,x70];%����:����ֵ,90���ŵ��Ͻ�,90���ŵ��½�,80���ŵ��Ͻ�,80���ŵ��½�,70���ŵ��Ͻ�,70���ŵ��½�
    end
    qwb{j}=qw;
end
save data10_deal